home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / cheb / init.c < prev    next >
Encoding:
C/C++ Source or Header  |  2002-04-18  |  2.5 KB  |  101 lines
  1. /* cheb/init.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. #include <config.h>
  21. #include <stdlib.h>
  22. #include <gsl/gsl_math.h>
  23. #include <gsl/gsl_errno.h>
  24. #include <gsl/gsl_chebyshev.h>
  25.  
  26. /*-*-*-*-*-*-*-*-*-*-*-* Allocators *-*-*-*-*-*-*-*-*-*-*-*/
  27.  
  28. gsl_cheb_series * 
  29. gsl_cheb_alloc(const size_t order)
  30. {
  31.   gsl_cheb_series * cs = (gsl_cheb_series *) malloc(sizeof(gsl_cheb_series));
  32.   
  33.   if(cs == 0) {
  34.     GSL_ERROR_VAL("failed to allocate gsl_cheb_series struct", GSL_ENOMEM, 0);
  35.   }
  36.   
  37.   cs->order    = order;
  38.   cs->order_sp = order;
  39.  
  40.   cs->c = (double *) malloc((order+1) * sizeof(double));
  41.  
  42.   if(cs->c == 0) {
  43.     GSL_ERROR_VAL("failed to allocate cheb coefficients", GSL_ENOMEM, 0);
  44.   }
  45.  
  46.   cs->f = (double *) malloc((order+1) * sizeof(double));
  47.  
  48.   if(cs->f == 0) {
  49.     GSL_ERROR_VAL("failed to allocate cheb function space", GSL_ENOMEM, 0);
  50.   }
  51.  
  52.   return cs;
  53. }
  54.  
  55.  
  56. void gsl_cheb_free(gsl_cheb_series * cs)
  57. {
  58.   free(cs->c);
  59.   free(cs);
  60. }
  61.  
  62. /*-*-*-*-*-*-*-*-*-*-*-* Initializer *-*-*-*-*-*-*-*-*-*-*-*/
  63.  
  64. int gsl_cheb_init(gsl_cheb_series * cs, const gsl_function *func,
  65.                   const double a, const double b)
  66. {
  67.   size_t k, j;
  68.  
  69.   if(a >= b) {
  70.     GSL_ERROR_VAL("null function interval [a,b]", GSL_EDOM, 0);
  71.   }
  72.   cs->a = a;
  73.   cs->b = b;
  74.   /* cs->err = 0.0; */
  75.  
  76.   { 
  77.     double bma = 0.5 * (cs->b - cs->a);
  78.     double bpa = 0.5 * (cs->b + cs->a);
  79.     double fac = 2.0/(cs->order +1.0);
  80.  
  81.     for(k = 0; k<=cs->order; k++) {
  82.       double y = cos(M_PI * (k+0.5)/(cs->order+1));
  83.       cs->f[k] = GSL_FN_EVAL(func, (y*bma + bpa));
  84.     }
  85.     
  86.     for(j = 0; j<=cs->order; j++) {
  87.       double sum = 0.0;
  88.       for(k = 0; k<=cs->order; k++) 
  89.         sum += cs->f[k]*cos(M_PI * j*(k+0.5)/(cs->order+1));
  90.       cs->c[j] = fac * sum;
  91.     }
  92.     
  93.   }
  94.   return GSL_SUCCESS;
  95. }
  96.  
  97.  
  98.  
  99.  
  100.  
  101.